{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "million-primary",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "import random\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "treated-universal",
   "metadata": {},
   "outputs": [],
   "source": [
    "tlr_miu = 0.001 # 求miu,s,b表的tolerance\n",
    "tlr_model = 0.0001 # 模型检验的tolerance\n",
    "precision = 2 # 参数保留小数点后的位数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cloudy-paradise",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 根据miu和s求b的函数\n",
    "def get_b(miu, s):\n",
    "    b = 0\n",
    "    for r in range(s):\n",
    "        b += (miu**r) * math.e**(-miu) / (math.factorial(r))\n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "worth-particular",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 单调递减函数二分法求变量值\n",
    "def binary_get(aim, func=get_b, s=0, lb=0, ub=1, tlr=tlr_miu, precision=3):\n",
    "    \"\"\"\n",
    "    aim:目标值y\n",
    "    func:y=f(x)的函数f\n",
    "    s:函数的额外参数\n",
    "    lb:x的下界\n",
    "    ub:x的上界\n",
    "    tlr:算法的精度，即y的精度\n",
    "    precision:返回的x保留几位小数\n",
    "    return:给定精度下目标值y在函数func下对应的x\n",
    "    \"\"\"\n",
    "    c = (lb + ub) / 2\n",
    "    valc = func(c, s)\n",
    "    if abs(valc - aim) < tlr:\n",
    "        return round(c, precision)\n",
    "    if aim < valc:\n",
    "        return binary_get(aim, func, s, c, ub, tlr, precision)\n",
    "    return binary_get(aim, func, s, lb, c, tlr, precision)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "03bde5c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_dist_per(x1, y1, x2, y2):\n",
    "    rady1 = math.radians(y1)\n",
    "    rady2 = math.radians(y2)\n",
    "    a = rady1 - rady2\n",
    "    b = math.radians(x1) - math.radians(x2)\n",
    "    s = 2*math.asin(math.sqrt(math.sin(a/2)**2 + math.cos(rady1)*math.cos(rady2)* math.sin(b/2)**2)) * 6378.004\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "preceding-bride",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获得距离矩阵\n",
    "def get_dist(group1, group2):\n",
    "    dist_mat = np.zeros([len(group1), len(group2)])\n",
    "    for i in range(len(group1)):\n",
    "        for j in range(len(group2)):\n",
    "            #dist_mat[i,j] = ((group1[i][0]-group2[j][0])**2 + (group1[i][1]-group2[j][1])**2)**0.5\n",
    "            dist_mat[i, j] = get_dist_per(group1[i][0], group1[i][1], group2[j][0], group2[j][1])\n",
    "    return dist_mat.round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "satisfactory-chance",
   "metadata": {},
   "source": [
    "### 参数设置"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "alpha-nomination",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_data = pd.read_excel(\"Test data.xlsx\", sheet_name='data')\n",
    "df_data_t = pd.read_excel('Test data.xlsx', sheet_name='data_t_jk')\n",
    "df_data_h = pd.read_excel('Test data.xlsx', sheet_name='data_h_jk')\n",
    "df_data_info_cust = pd.read_excel('Test data.xlsx', sheet_name='info_cust')\n",
    "with open('data1.txt', \"r\", encoding='utf-8') as fl:\n",
    "    C_num = int(fl.readline().split()[0])\n",
    "    J_num = int(fl.readline().split()[0])\n",
    "    I_num = int(fl.readline().split()[0])\n",
    "    M_num = int(fl.readline().split()[0])\n",
    "    virtual_num = int(fl.readline().split()[0])\n",
    "    K_num = int(fl.readline().split()[0])\n",
    "    per_cost = list(df_data.iloc[:, 1].values)\n",
    "    # per_cost = list(map(int, fl.readline().split()[0].split(\",\")))\n",
    "    h = df_data_h.values[:, 1:].astype(\"int32\")\n",
    "    \"\"\"temp = []\n",
    "    for line in fl.readline().split()[0].split(\";\"):\n",
    "        temp.append(list(map(int, line.split(\",\"))))\n",
    "    h = np.array(temp)\"\"\"\n",
    "    t = df_data_t.values[:, 1:].astype(\"int32\")\n",
    "    \"\"\"temp = []\n",
    "    for line in fl.readline().split()[0].split(\";\"):\n",
    "        temp.append(list(map(int, line.split(\",\"))))\n",
    "    t = np.array(temp)\"\"\"\n",
    "    v_1 = int(fl.readline().split()[0])\n",
    "    v_2 = list(df_data.iloc[:, 2].values.astype('int'))\n",
    "    #v_2 = list(map(int, fl.readline().split()[0].split(\",\")))\n",
    "    max_s = int(fl.readline().split()[0])\n",
    "    virtual_per_cost = list(df_data.iloc[:, 1].values.astype('int'))\n",
    "    #virtual_per_cost = list(map(int, fl.readline().split()[0].split(\",\")))\n",
    "\n",
    "with open('info_center.txt', \"r\", encoding='utf-8') as fl:\n",
    "    c_coords = []\n",
    "    for line in fl.readlines()[1:]:\n",
    "        c_coords.append(tuple(map(float, line.split()[1:])))\n",
    "\n",
    "with open('info_depot.txt', \"r\", encoding='utf-8') as fl:\n",
    "    j_coords = []\n",
    "    f = []\n",
    "    for line in fl.readlines()[1:]:\n",
    "        j_coords.append(tuple(map(float, line.split()[1:3])))\n",
    "        f.append(int(line.split()[-1]))\n",
    "\n",
    "with open('info_express.txt', \"r\", encoding='utf-8') as fl:\n",
    "    m_coords = []\n",
    "    t_res = []\n",
    "    for line in fl.readlines()[1:]:\n",
    "        m_coords.append(tuple(map(float, line.split()[1:3])))\n",
    "        t_res.append(int(line.split()[-1]))\n",
    "    for coord in j_coords:\n",
    "        m_coords.append(coord)\n",
    "        t_res.append(0)\n",
    "\n",
    "with open('info_cust1.txt', \"r\", encoding='utf-8') as fl:\n",
    "    i_coords = []\n",
    "    D = []\n",
    "    alpha = []\n",
    "    T = []\n",
    "    for line in fl.readlines()[1:]:\n",
    "        temp = line.split()\n",
    "        i_coords.append(tuple(map(float, temp[1:3])))\n",
    "        \"\"\"D.append(list(map(float, temp[3].split(\",\"))))\n",
    "        alpha.append(list(map(float, temp[4].split(\",\"))))\n",
    "        T.append(list(map(float, temp[5].split(\",\"))))\"\"\"\n",
    "    \"\"\"D = np.array(D)\n",
    "    alpha = np.array(alpha)\n",
    "    T = np.array(T)\"\"\"\n",
    "    D = np.zeros((I_num, K_num))\n",
    "    alpha = np.zeros((I_num, K_num))\n",
    "    T = np.zeros((I_num, K_num))\n",
    "    for _, cust_n, prod_n, d_n, alpha_n, t_n in df_data_info_cust.itertuples(index=False):\n",
    "        D[cust_n,prod_n] = d_n\n",
    "        alpha[cust_n,prod_n] = alpha_n\n",
    "        T[cust_n,prod_n] = t_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "buried-redhead",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.88, 0.75, 0.82, 0.85, 0.9 , 0.78, 0.91, 0.88, 0.81, 0.8 , 0.78,\n",
       "        0.86, 0.85, 0.84, 0.84, 0.85, 0.8 , 0.8 , 0.9 , 0.82],\n",
       "       [0.82, 0.82, 0.76, 0.76, 0.78, 0.9 , 0.82, 0.84, 0.86, 0.85, 0.84,\n",
       "        0.88, 0.8 , 0.78, 0.9 , 0.89, 0.88, 0.9 , 0.86, 0.78],\n",
       "       [0.76, 0.78, 0.88, 0.82, 0.84, 0.86, 0.84, 0.86, 0.88, 0.9 , 0.92,\n",
       "        0.88, 0.86, 0.82, 0.82, 0.9 , 0.82, 0.82, 0.84, 0.86],\n",
       "       [0.84, 0.86, 0.85, 0.88, 0.89, 0.9 , 0.9 , 0.85, 0.88, 0.84, 0.88,\n",
       "        0.88, 0.88, 0.88, 0.88, 0.88, 0.86, 0.88, 0.88, 0.9 ],\n",
       "       [0.88, 0.88, 0.88, 0.88, 0.86, 0.75, 0.88, 0.88, 0.88, 0.88, 0.9 ,\n",
       "        0.88, 0.88, 0.86, 0.88, 0.88, 0.88, 0.86, 0.88, 0.88],\n",
       "       [0.88, 0.88, 0.88, 0.9 , 0.88, 0.88, 0.84, 0.88, 0.88, 0.88, 0.88,\n",
       "        0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88],\n",
       "       [0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88,\n",
       "        0.84, 0.88, 0.88, 0.88, 0.86, 0.88, 0.88, 0.88, 0.88],\n",
       "       [0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.83, 0.88, 0.88, 0.88, 0.88,\n",
       "        0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88],\n",
       "       [0.9 , 0.88, 0.88, 0.88, 0.88, 0.88, 0.85, 0.88, 0.85, 0.88, 0.88,\n",
       "        0.76, 0.88, 0.78, 0.88, 0.8 , 0.82, 0.88, 0.84, 0.88],\n",
       "       [0.84, 0.88, 0.76, 0.88, 0.78, 0.88, 0.88, 0.88, 0.84, 0.82, 0.88,\n",
       "        0.88, 0.86, 0.88, 0.74, 0.78, 0.88, 0.84, 0.78, 0.88]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "welcome-elements",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"C_num = 1 # 中心库数量\n",
    "J_num = 3 # 本地库备选点数量\n",
    "I_num = 10 # 客户数量\n",
    "M_num = 8 # 快运站集合（包含虚拟点）\n",
    "virtual_num = J_num # 虚拟节点的数量，位置和本地库本选点一样\n",
    "K_num = 3 # 备件数量\n",
    "\n",
    "per_cost = [random.randint(1,6)*0.1 for _ in range(K_num)]\"\"\"\n",
    "\n",
    "C = {i for i in range(C_num)}\n",
    "I = {i for i in range(I_num)}\n",
    "J = {i for i in range(J_num)}\n",
    "M = {i for i in range(M_num)}\n",
    "K = {i for i in range(K_num)}\n",
    "\n",
    "\"\"\"c_coords = [(random.randint(0,100), random.randint(0,100))]\n",
    "i_coords = [(random.randint(0,100), random.randint(0,100)) for _ in range(len(I))]\n",
    "j_coords = [(random.randint(0,100), random.randint(0,100)) for _ in range(len(J))]\n",
    "m_coords = [(random.randint(0,100), random.randint(0,100)) for _ in range(len(M)-virtual_num)] + j_coords\n",
    "\"\"\"\n",
    "dist_ci = get_dist(c_coords, i_coords)\n",
    "dist_jm = get_dist(j_coords, m_coords)\n",
    "dist_ji = get_dist(j_coords, i_coords)\n",
    "\n",
    "\"\"\"f = []\n",
    "for j in range(len(J)):\n",
    "    f.append(random.randint(500,1500))\"\"\"\n",
    "    \n",
    "c_4 = np.zeros([len(I), len(M), len(J), len(K)])\n",
    "for i in range(len(c_4)):\n",
    "    for m in range(len(c_4[0])):\n",
    "        for j in range(len(c_4[0,0])):\n",
    "            for k in range(len(c_4[0,0,0])):\n",
    "                if m < M_num - virtual_num:\n",
    "                    c_4[i,m,j,k] = (dist_jm[j,m] + dist_ji[j,i]) * per_cost[k]\n",
    "                else:\n",
    "                    c_4[i,m,j,k] = (dist_jm[j,m] + dist_ji[j,i]) * virtual_per_cost[k]\n",
    "\n",
    "    \n",
    "c_3 = np.zeros([len(C), len(I), len(K)])\n",
    "for c in C:\n",
    "    for i in I:\n",
    "        for k in K:\n",
    "            c_3[c,i,k] = dist_ci[c,i] * virtual_per_cost[k]\n",
    "\n",
    "d_jm = dist_jm\n",
    "d_ji = dist_ji\n",
    "d_ci = dist_ci\n",
    "\n",
    "\"h = np.random.randint(10,20,size=(len(J), len(K)))\"\n",
    "\n",
    "\"D = np.random.randint(0,2,size=(len(I),len(K)))\"\n",
    "\n",
    "\"t = np.random.randint(1,3, size=(len(J), len(K)))\"\n",
    "\n",
    "\"\"\"v_1 = 1 # 本地库-快运点的行驶速度\n",
    "v_2 = [random.randint(10,20)*0.1 for _ in K]\"\"\"\n",
    "tao_jm = d_jm/v_1\n",
    "tao_jik = np.zeros((len(J), len(I), len(K)))\n",
    "for k in K:\n",
    "    tao_jik[:,:,k] = (d_ji/v_2[k]).round(2)\n",
    "\n",
    "\"alpha = np.random.random((len(I), len(K))).round(2)*0.5\"\n",
    "\n",
    "\"t_res = [random.randint(1,10)*0.1+1 for _ in M] # 快运点m的服务相应时间\"\n",
    "\n",
    "\"T = np.ones((len(I), len(K))) * 80\"\n",
    "\n",
    "\"max_s = 5\"\n",
    "L = dict() # 可用库存水平集合\n",
    "for j in J:\n",
    "    for k in K:\n",
    "        L[j,k] = set(i+1 for i in range(max_s))\n",
    "\n",
    "omega = T # 客户i对产品k要求的运输时间限制\n",
    "\n",
    "miu_lb = 0\n",
    "miu_ub = np.sum(D)*10\n",
    "\n",
    "J_ik = dict() # 能服务客户i的所有本地库备选点集合\n",
    "I_jk = dict() # 本地库j能服务的所有客户集合\n",
    "for i in I:\n",
    "    for j in J:\n",
    "        for k in K:\n",
    "            if (i,k) not in J_ik.keys():\n",
    "                J_ik[i,k] = set()\n",
    "            if (j,k) not in I_jk.keys():\n",
    "                I_jk[j,k] = set()\n",
    "            if tao_jik[j,i,k] < omega[i,k]:\n",
    "                J_ik[i,k].add(j)\n",
    "                I_jk[j,k].add(i)\n",
    "\n",
    "ijk_b = dict()\n",
    "bjk_i = dict()\n",
    "for i in I:\n",
    "    for j in J:\n",
    "        for k in K:\n",
    "            ijk_b[i,j,k] = alpha[i,k]\n",
    "for i in I:\n",
    "    for j in J:\n",
    "        for k in K:\n",
    "            b = ijk_b[i,j,k]\n",
    "            if (b,j,k) not in bjk_i.keys():\n",
    "                bjk_i[b,j,k] = i\n",
    "                \n",
    "\n",
    "                \n",
    "MC =1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "advised-illinois",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 给定一个b,i,j,k,s，求对应的miu值\n",
    "miu = dict()\n",
    "for i in I:\n",
    "    for j in J:\n",
    "        for k in K:\n",
    "            b = ijk_b[i,j,k]\n",
    "            for s in L[j,k]:\n",
    "                miu[i,j,k,s] = binary_get(b,s=s,lb=miu_lb, ub=miu_ub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "tough-royalty",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = -1\n",
    "j = -1\n",
    "k = -1\n",
    "m = -1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "taken-reaction",
   "metadata": {},
   "source": [
    "### 建模"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "played-liabilities",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file C:\\Users\\CYK\\gurobi.lic\n"
     ]
    }
   ],
   "source": [
    "m = gp.Model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "quantitative-allergy",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = m.addVars(len(J), vtype=GRB.BINARY)\n",
    "x = m.addVars(len(I), len(M), len(J), len(K), vtype=GRB.BINARY)\n",
    "v = m.addVars(len(J), max_s+1, len(K), vtype=GRB.BINARY)\n",
    "w = m.addVars(len(J), len(I), len(K),vtype=GRB.BINARY)\n",
    "z = m.addVars(len(I), len(M), len(J), len(I), len(K), vtype=GRB.BINARY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f9a7a93a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# z值确定式\n",
    "m.addConstrs(z[i,m,j,n,k] - x[i,m,j,k] - w[j,n,k] >= -1 for i in I for m in M for j in J for n in I for k in K)\n",
    "m.addConstrs(z[i,m,j,n,k] - x[i,m,j,k] - w[j,n,k] <= 1 for i in I for m in M for j in J for n in I for k in K)\n",
    "m.addConstrs(z[i,m,j,n,k] - x[i,m,j,k] <= 0 for i in I for m in M for j in J for n in I for k in K)\n",
    "m.addConstrs(z[i,m,j,n,k] - w[j,n,k] <= 0 for i in I for m in M for j in J for n in I for k in K)\n",
    "print(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "presidential-pride",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 当本地库开放时，才可以进行分配\n",
    "m.addConstrs(x[i,m,j,k] - y[j] <= 0 for i in I for k in K for m in M for j in J_ik[i,k])\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "tutorial-expert",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 每个客户需求必须且仅能分配给一个本地库和一个快运点\n",
    "m.addConstrs(gp.quicksum(gp.quicksum(x[i,m,j,k] for m in M) for j in J_ik[i,k]) == 1 for i in I for k in K)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "indonesian-newspaper",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 运输时间约束\n",
    "m.addConstrs(gp.quicksum(tao_jm[j,m]*x[i,m,j,k] for j in J_ik[i,k]) + \\\n",
    "            gp.quicksum(tao_jik[j,i,k]*x[i,m,j,k] for j in J_ik[i,k]) <= \\\n",
    "             T[i,k] - t_res[m] for i in I for k in K for m in M)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "sound-liberty",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 服务水平约束\n",
    "m.addConstrs(gp.quicksum(ijk_b[n,j,k]*w[j,n,k] for n in I) - gp.quicksum(alpha[i,k]*x[i,m,j,k] for m in M) >= 0 for i in I for k in K for j in J_ik[i,k])\n",
    "m.addConstrs(gp.quicksum(w[j,n,k] for n in I) <= 1 for j in J for k in K)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "concrete-socket",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 库存水平约束\n",
    "m.addConstrs(gp.quicksum(v[j,s,k] for s in L[j,k]) <= 1 for j in J for k in K)\n",
    "m.addConstrs(t[j,k] * gp.quicksum(gp.quicksum(D[i,k]*x[i,m,j,k] for i in I_jk[j,k]) for m in M) - \\\n",
    "            gp.quicksum(miu[n,j,k,s]*v[j,s,k] for s in L[j,k]) + MC*w[j,n,k] <= MC for j in J for k in K for n in I)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4d04a5df",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)\n",
      "Thread count: 4 physical cores, 4 logical processors, using up to 4 threads\n",
      "Optimize a model with 367472 rows, 101463 columns and 947820 nonzeros\n",
      "Model fingerprint: 0x4e733563\n",
      "Variable types: 0 continuous, 101463 integer (101463 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [8e-02, 1e+03]\n",
      "  Objective range  [6e+01, 1e+05]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 1e+03]\n",
      "Presolve time: 0.04s\n",
      "\n",
      "Explored 0 nodes (0 simplex iterations) in 0.11 seconds\n",
      "Thread count was 1 (of 4 available processors)\n",
      "\n",
      "Solution count 0\n",
      "\n",
      "Model is infeasible\n",
      "Best objective -, best bound -, gap -\n"
     ]
    }
   ],
   "source": [
    "# 原模型\n",
    "lin1 = gp.quicksum(f[j]*y[j] for j in J)\n",
    "lin2 = gp.quicksum(gp.quicksum(gp.quicksum(h[j,k]*s*v[j,s,k] for s in L[j,k]) for k in K) for j in J)\n",
    "lin3 = gp.quicksum(gp.quicksum(gp.quicksum(gp.quicksum(gp.quicksum((c_4[i,m,j,k])*z[i,m,j,n,k]*D[i,k] for i in I_jk[j,k]) for m in M) for k in K) for j in J) for n in I)\n",
    "m.setObjective(lin1 + lin2 + lin3, sense=GRB.MINIMIZE)\n",
    "m.update()\n",
    "m.optimize()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f87d0b01",
   "metadata": {},
   "source": [
    "# 不考虑期望成本\n",
    "lin1 = gp.quicksum(f[j]*y[j] for j in J)\n",
    "lin2 = gp.quicksum(gp.quicksum(gp.quicksum(h[j,k]*s*v[j,s,k] for s in L[j,k]) for k in K) for j in J)\n",
    "lin3 = gp.quicksum(gp.quicksum(gp.quicksum(gp.quicksum(c_4[i,m,j,k]*x[i,m,j,k]*D[i,k] for i in I_jk[j,k]) for m in M) for k in K) for j in J)\n",
    "m.setObjective(lin1 + lin2 + lin3, sense=GRB.MINIMIZE)\n",
    "m.update()\n",
    "m.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "geographic-anatomy",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<gurobi.Model MIP instance Unnamed: 367472 constrs, 101463 vars, No parameter changes>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fca19d8",
   "metadata": {},
   "source": [
    "### 选址成本，库存成本，运输成本"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "23c71667",
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "Unable to retrieve attribute 'x'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-22-fd3ee9cdb355>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# 原模型\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mCost_location\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mCost_inventory\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0ms\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mL\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mK\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mCost_trans\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mc_4\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mz\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mD\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mI_jk\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mm\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mM\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mK\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mn\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mI\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mCost_location\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mCost_inventory\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mCost_trans\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-22-fd3ee9cdb355>\u001b[0m in \u001b[0;36m<genexpr>\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# 原模型\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mCost_location\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mCost_inventory\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0ms\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mL\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mK\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mCost_trans\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mc_4\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mz\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mD\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mI_jk\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mm\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mM\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mK\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mn\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mI\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mCost_location\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mCost_inventory\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mCost_trans\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32msrc\\gurobipy\\var.pxi\u001b[0m in \u001b[0;36mgurobipy.Var.__getattr__\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32msrc\\gurobipy\\var.pxi\u001b[0m in \u001b[0;36mgurobipy.Var.getAttr\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32msrc\\gurobipy\\attrutil.pxi\u001b[0m in \u001b[0;36mgurobipy.__getattr\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: Unable to retrieve attribute 'x'"
     ]
    }
   ],
   "source": [
    "# 原模型\n",
    "Cost_location = sum(f[j]*y[j].x for j in J)\n",
    "Cost_inventory = sum(sum(sum(h[j,k]*s*v[j,s,k].x for s in L[j,k]) for k in K) for j in J)\n",
    "Cost_trans = sum(sum(sum(sum(sum((c_4[i,m,j,k])*z[i,m,j,n,k].x*D[i,k] for i in I_jk[j,k]) for m in M) for k in K) for j in J) for n in I)\n",
    "print(Cost_location,Cost_inventory,Cost_trans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "corporate-dynamics",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Cost_location + Cost_inventory + Cost_trans"
   ]
  },
  {
   "cell_type": "raw",
   "id": "74bbceed",
   "metadata": {},
   "source": [
    "# 不考虑期望成本\n",
    "Cost_location = sum(f[j]*y[j].x for j in J)\n",
    "Cost_inventory = sum(sum(sum(h[j,k]*s*v[j,s,k].x for s in L[j,k]) for k in K) for j in J)\n",
    "Cost_trans = sum(sum(sum(sum(c_4[i,m,j,k]*x[i,m,j,k].x*D[i,k] for i in I_jk[j,k]) for m in M) for k in K) for j in J)\n",
    "print(Cost_location,Cost_inventory,Cost_trans)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "excited-functionality",
   "metadata": {},
   "source": [
    "### 检验\n",
    "##### 1. 每个客户必须且仅能分配给一个快递点和一个本地库\n",
    "##### 2. 保证被本地库j服务的所有客户都必须满足客户的服务水平\n",
    "##### 3. 库存水平约束"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "three-percentage",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1\n",
    "covers = dict() # 备选点j cover的（客户j, 备件k）\n",
    "routes_fori = dict() # （客户i，备件k）经过的（快运点m，本地库j）\n",
    "for i in I:\n",
    "    for m in M:\n",
    "        for j in J:\n",
    "            for k in K:\n",
    "                if abs(x[i,m,j,k].x - 1) < tlr_model:\n",
    "                    if j not in covers:\n",
    "                        covers[j] = set()\n",
    "                    covers[j].add((i,k))\n",
    "                    if (i,k) in routes_fori:\n",
    "                        raise \"错误\"\n",
    "                    routes_fori[i,k] = (m,j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "instructional-theorem",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc10fe9b",
   "metadata": {},
   "source": [
    "### 查询结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "british-failing",
   "metadata": {},
   "outputs": [],
   "source": [
    "tlr = 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "specific-elder",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 开放方案\n",
    "for j in J:\n",
    "    if abs(y[j].x - 1) < tlr:\n",
    "        print(j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vulnerable-journalism",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 顾客-产品分配方案\n",
    "ll = []\n",
    "have_error = {}\n",
    "for i in I:\n",
    "    for k in K:\n",
    "        for m in M:\n",
    "            for j in J:\n",
    "                if abs(x[i,m,j,k].x - 1) < tlr:\n",
    "                    if x[i,m,j,k].x != 1:\n",
    "                        have_error[i,m,j,k] = x[i,m,j,k].x\n",
    "                    if D[i,k] == 0:\n",
    "                        continue\n",
    "                    if m >= M_num - virtual_num:\n",
    "                        ll.append(-1)\n",
    "                        print('({0}, {1})分配的仓库是{2},运输方式是直运'.format(i,k,j))\n",
    "                    else:\n",
    "                        ll.append(m)\n",
    "                        print('({0}, {1})分配的仓库是{2},运输方式是{3}'.format(i,k,j,m))\n",
    "print(ll)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "infrared-cameroon",
   "metadata": {},
   "outputs": [],
   "source": [
    "error = {}\n",
    "for i in I:\n",
    "    for m in M:\n",
    "        for j in J:\n",
    "            for n in I:\n",
    "                for k in K:\n",
    "                    if abs(z[i,m,j,n,k].x -1) < tlr and z[i,m,j,n,k].x != 1:\n",
    "                        error[i,m,j,n,k] = z[i,m,j,n,k].x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hydraulic-likelihood",
   "metadata": {},
   "outputs": [],
   "source": [
    "error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "absolute-nebraska",
   "metadata": {},
   "outputs": [],
   "source": [
    "j=1\n",
    "m=1\n",
    "k=1\n",
    "n=9\n",
    "i=6\n",
    "error_t = (c_4[i,m,j,k]*ijk_b[n,j,k] + c_3[0,i,k]*(1-ijk_b[n,j,k]))*(1-z[i,m,j,n,k].x)*D[i,k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "periodic-friday",
   "metadata": {},
   "outputs": [],
   "source": [
    "Cost_trans+error_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "broad-desperate",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 库存水平\n",
    "for j in J:\n",
    "    for s in range(max_s+1):\n",
    "        for k in K:\n",
    "            if abs(v[j,s,k].x - 1) < tlr:\n",
    "                print(\"仓库{0}对产品{1}的库存水平为：{2}\".format(j,k,s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "british-explorer",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argmax(alpha[:, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "minor-constraint",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 9\n",
    "j = 1\n",
    "k = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "isolated-representation",
   "metadata": {},
   "outputs": [],
   "source": [
    "t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "later-enhancement",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "miu[i,j,k,7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "gentle-click",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "comic-ensemble",
   "metadata": {},
   "outputs": [],
   "source": [
    "miu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fourth-packaging",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(D[:, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "excessive-thousand",
   "metadata": {},
   "outputs": [],
   "source": [
    "D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d972ba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "tt = 0.000001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "957b0540",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for j in range(len(J)):\n",
    "    if abs(y[j].x - 1)< tt:\n",
    "        print(j)\n",
    "for i in range(len(I)):\n",
    "    for m in range(len(M)):\n",
    "        for j in range(len(J)):\n",
    "            for k in range(len(K)):\n",
    "                if abs(x[i,m,j,k].x - 1) < tt:\n",
    "                    print(i,m,j,k)\n",
    "for j in range(len(J)):\n",
    "    for s in range(max_s+1):\n",
    "        for k in range(len(K)):\n",
    "            if abs(v[j,s,k].x - 1) < tt:\n",
    "                print(j,s,k)\n",
    "for j in range(len(J)):\n",
    "    for i in range(len(I)):\n",
    "        for k in range(len(K)):\n",
    "            if abs(w[j,i,k].x -1) < tt:\n",
    "                print(j,i,k)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "3b9a8c2a",
   "metadata": {},
   "source": [
    "for j in range(len(J)):\n",
    "    if abs(y[j].x - 1) > tt and abs(y[j].x) > tt:\n",
    "        print(j)\n",
    "for i in range(len(I)):\n",
    "    for m in range(len(M)):\n",
    "        for j in range(len(J)):\n",
    "            for k in range(len(K)):\n",
    "                if abs(x[i,m,j,k].x - 1) > tt and abs(x[i,m,j,k].x) > tt:\n",
    "                    print(i,m,j,k)\n",
    "for j in range(len(J)):\n",
    "    for s in range(max_s+1):\n",
    "        for k in range(len(K)):\n",
    "            if abs(v[j,s,k].x - 1) > tt and abs(v[j,s,k].x) > tt:\n",
    "                print(j,s,k)\n",
    "for j in range(len(J)):\n",
    "    for i in range(len(I)):\n",
    "        for k in range(len(K)):\n",
    "            if abs(w[j,i,k].x -1) > tt and abs(w[j,i,k].x) > tt:\n",
    "                print(j,i,k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "listed-chaos",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "moved-dispute",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "casual-mathematics",
   "metadata": {},
   "outputs": [],
   "source": [
    "fl = open('xx.txt', 'w')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "classical-northern",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(fl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "complete-locking",
   "metadata": {},
   "outputs": [],
   "source": [
    "fl.__class__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "excessive-delaware",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
